#ifndef PBVP_BSPLINE_INTERPOL_H
#define PBVP_BSPLINE_INTERPOL_H

/************************************************************************/
/* BSpline Interpolation of volume functions
   Based on Xu Guoliang's code
   
   Notice the interpolation and derivatives are not accurate near the 
   boundary of the volume.
   */
/************************************************************************/

namespace PBVP {
	class BSplineInterpol {
	public:
		/// Convert a 1D function to interpolating B-SPline coefficients
		/** 
		 * @param s input samples --> output coefficients 
		 * @param size	number of data(coefficients)
		 * @param epsilon error tolerance
		 */
		static void interpolCoeff1D(float *s, int size, float epsilon = DEFAULT_BSI_EPSILON);

		/// Convert a 2D function to interpolating B-Spline coefficients
		/** 
		* @param s input samples --> output coefficients 
		* @param nx	number of samples(coefficients) in the x direction
		* @param ny number of samples(coefficients) in the y direction
		* @param epsilon error tolerance
		*/
		static void interpolCoeff2D(float *s, int nx, int ny, float epsilon = DEFAULT_BSI_EPSILON);


		/// Convert a 3D function to interpolating B-Spline coefficients
		/** 
		* @param s input samples --> output coefficients 
		* @param nx	number of samples(coefficients) in the x direction
		* @param ny number of samples(coefficients) in the y direction
		* @param nz number of samples(coefficients) in the z direction
		* @param epsilon error tolerance
		*/
		static void interpolCoeff3D(float *s, int nx, int ny, int nz, float epsilon = DEFAULT_BSI_EPSILON);

		/// Compute degree 1 and 2 partial derivatives of a B-Spline interpolated function
		/** 
		 * @param c The spline coefficients
		 * @param dx spacing in the x direction
		 * @param dy spacing in the y direction
		 * @param dz spacing in the z direction
		 * @param nx number of samples(coefficients) in the x direction
		 * @param ny number of samples(coefficients) in the y direction
		 * @param nz number of samples(coefficients) in the z direction
		 * @param partials A 10-element array containing returned 1st and 2nd order partial derivatives \
		 *				   The 0th element is the func value, elements 1-3 are 1st derivatives \
		 *				   elements 4 - 9 are 2nd derivatives
		 */
		static void partialDeriv12(float *c, float dx, float dy, float dz,
							int nx, int ny, int nz,
							int u, int v, int w, float *partials);

		///
		static void computeTensor();

	private:
		static const float DEFAULT_BSI_EPSILON;
		static float TensorF[27], TensorFx[27], TensorFy[27], TensorFz[27], TensorFxx[27],
					 TensorFxy[27], TensorFxz[27], TensorFyy[27], TensorFyz[27], TensorFzz[27];
	
		///
		static void Tensor_333(float *xx, float *yy, float *zz, float *result);

		/// 
		static void get27Coefficients(float *c, int nx, int ny, int nz, int u, int v, int w, float *c27);

		/// 
		static float getOneCoefficient(float *c, int nx, int ny, int nz, int u, int v, int w); 

	};

}

#endif